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ABSTRACT 


THE  BOUNDARY  INTEGRAL  EQUATION  METHOD  USING  VARIOUS 
APPROXIMATION  TECHNIQUES  FOR  PROBLEMS  GOVERNED  BY 

LAPLACE'S  EQUATION 

The  Boundary  Integral  Equation  (BIE)  method  is  applied  to  boun- 
dary value  problems  governed  by  Laplace's  equation*  In  addition  to 
the  piecewise  constant  approximation  for  boundary  functions,  two  other 
numerical  approximations,  quadratic  shape  function  and  cubic  spline 
function  representations,  are  adopted.  Comparisons  are  discussed. 
Topics  for  further  research  are  indicated. 
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CHAPTER  ONE 
INTRODUCTION 

Probably  the  most  familiar  partial  differential  equation  of  the 
elliptic  type  is  Laplace's  equation: 

V2U  = 0 

2 2 2 
2 ri  a O 

where  v = =— — + — - + — ~ for  Cartesian  coordinates. 

2 b Lt 

bx  by  b z 

It  plays  an  important  role  in  the  various  branches  of  the  engineering 
and  applied  mathematics.  It  is  satisfied,  for  example,  if  U is 

(a)  the  temperature  distribution  under  steady  conditions  in  a region 
of  constant  conductivity; 

(b)  the  warping  function  in  the  torsion  of  cylindrical  bars; 

(c)  either  the  stream  function  or  the  velocity  potential  in  the  irrota- 
tional  flow  of  an  incompressible  inviscid  fluid; 

(d)  any  distribution  of  electric  potential  in  a region  of  constant  re- 
sistivity; 

(e)  the  distribution  of  magnetic  potential  in  regions  of  constant  per- 
meability; 

(f)  either  the  real  or  imaginary  part  of  any  analytic  function  of  a 
complex  variable; 

Therefore,  the  treatment  of  solutions  to  boundary  value  problems  for 
this  differential  equation  is  interesting  and  important. 

The  methods  of  solution  for  such  problems  include  analytical  methods, 
graphical  methods,  numerical  methods  and  analogical  methods.  For 
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the  analytical  approach,  the  method  of  separation  of  variables  is  used 
widely,  where  the  general  pattern  yields  different  forms  of  solutions 
depending  upon  the  nature  of  the  problem  and  the  particular  coordinate 
systems  involved.  However,  it  is  not  always  possible  to  separate  vari- 
ables [l]*.  As  to  another  analytical  approach,  the  similarity  method, 
the  similarity  transformation  [2]  can  be  used  via  group  theory  such 
that  Laplace's  equation  will  be  reduced  to  an  ordinary  differential 
equation  easily.  However,  the  capability  of  this  transformation  to  pro- 
perly carry  over  the  boundary  conditions  to  be  compatible  with  the 
transformed  equation  is  not  always  guaranteed.  In  fact,  even  though 
there  are  some  other  analytical  approaches,  e.g.,  the  conformal  map- 
ping method  [l],  it  turns  out  that  the  analytical  approach  is  recom- 
mended only  for  dealing  with  problems  which  are  simple  in  geometry 
and  boundary  conditions.  When  the  geometry  or  boundary  conditions 
become  complicated,  the  analytical  approach  will  be  too  involved  to  be 
practical. 

For  the  graphical  method  [l],  even  with  complex  geometry,  the 
problem  will  be  easy  to  handle  if  the  boundary  conditions  are  uniform 
and  normal  derivatives  of  the  potential  function  itself  across  the  sur- 
face are  zero.  However,  it  will  be  awkward  if  the  boundary  conditions 
involved  are  not  uniform  or  there  are  non-zero  normal  derivatives  of 
the  potential  function  across  the  surface.  And  even  with  simple 
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problems,  the  error  of  results  can  be  very  large  without  very  great  care 
in  drawing. 

The  analogical  method  [3]  can  be  applied  to  complicated  problems, 
but  sometimes  the  cost  of  construction  of  an  analogue,  obtaining  the 
solutions  in  the  analogous  system,  and  translating  them  back  into  the 
original  is  tedious  and  not  straightforward. 

By  considering  all  fhe  advantages  and  disadvantages  of  the  differ- 
ent approaches  mentioned  above,  it  is  found  by  comparison  that  the 
numerical  technique  is  generally  the  most  desirable  solution  method, 
since  it  is  flexible  and  is  even  applicable  to  systems  with  variable  phy- 
sical properties  and  nonuniform  boundary  conditions.  More  important, 
in  recent  years  high-speed  and  large -capacity  digital  computers  have 
developed  rapidly  and  have  become  especially  suitable  for  complicated 
and  tedious  numerical  manipulations.  Therefore  the  importance  of  the 
numerical  method  has  grown  greatly  and  will  undoubtedly  continue  to 
grow  at  an  increasing  rate.  Moreover,  although  the  solution  of  prob- 
lems using  numerical  methods  does  not  provide  general  solutions  as 
do  the  analytical  methods,  there  are  many  situations  for  which  numeri- 
cal solutions  provide  the  only  approach  which  can  be  used.  In  fact,  it 
is  often  found  that  when  a general  solution  is  known,  it  proves  to  be 
very  difficult  and  tedious  to  translate  into  particular  results  for  a 

I 

particular  problem;  therefore,  we  may  say  that  not  only  are  numerical 
methods  essential  in  problems  which  will  not  yield  to  any  other  method 

t 

of  solution,  but  that  also  they  are  often  the  best  for  obtaining  a 


i1 

ij 


-4- 


I 


particular  solution  even  when  a general  solution  is  available  by  analyti- 
cal methods.  Moreover,  it  is  usually  the  particular  solution  which  the 
engineer  and  the  physicist  are  interested  in  obtaining. 

In  the  past,  various  numerical  techniques  have  been  applied  to 
solve  Laplacian  nroblems.  Among  them,  finite  difference  methods 
(FDM)  and  finite  element  methods  (FEM)  are  popular,  the  practical 
difficulties  associated  with  them  being  also  well-known  [4]  [5],  How- 
ever, by  noticing  that  Laplace's  equation  is  an  elliptic  type  partial 
differential  equation  and  therefore  Laplacian  problems  are  boundary 
value  problems,  a basic  question  will  be  raised:  Is  it  possible  to  find 
a numerical  technique  which  deals  as  much  as  possible  only  with  the 
boundary  data  instead  of  referring  at  the  outset  to  the  whole  domain, 
boundary  and  interior,  as  have  been  done  in  both  FDM  and  FEM  ? 

The  answer  is  affirmative.  Specifically,  the  boundary  integral  equa- 
tion method  (BIE)  [6]  is  appropriate  for  this  situation.  The  major 
feature  of  BIE  is  that  unknown  boundary  data  may  be  obtained  by  refer- 
ence only  to  the  boundary  of  the  domain  of  interest  without  involving 
the  interior.  In  essence,  a given  problem  is  solved  on  the  boundary 
and  the  desired  potential  is  evaluated  at  a desired  interior  point  in 
terms  of  the  boundary  data.  Obviously,  the  dimensions  of  the  problem 
have  been  reduced  by  one.  The  great  advantage  to  the  analyst  comes 
from  simplification  of  the  mathematical  modelling.  Thereafter,  in 
this  thesis,  the  BIE  method  will  be  applied  to  problems  for  Laplace's 
equation.  Different  approximation  techniques  used  to  carry  out  the 
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numerical  analysis,  having  different  features  and  advantages,  will  also 
be  discussed. 

In  Chapter  2 we  formulate  boundary  value  problems  for  Laplace's 
equation,  via  the  BIE,  for  two  and  three  dimensional  problems.  Also 
we  indicate  a basic  numerical  procedure  used  with  the  BIE. 

In  Chapter  3 we  discuss  improved  approximation  methods,  1.  e. 
the  quadratic  shape  function  approximation  for  two  and  three  dimen- 
sional problems  and  the  spline  function  approximation  for  two  dimen- 
sional problems.  All  procedures  necessary  to  integrate  the  various 
kernel  functions  to  form  the  elements  of  the  matrix  of  coefficients  in 
the  resultant  system  of  algebraic  equations  are  presented  in  three 
Appendices.  These  procedures,  though  given  in  the  Appendices  to  not 
interrupt  the  flow  of  the  formulative  ideas,  are  not  merely  detail.  In- 
deed these  procedures  form  the  basis  for  actual  numerical  solutions  as 
illustrated  in  section  3 of  Chapter  3. 

In  Chapter  4 we  offer  some  discussion  and  conclusions  about  the 
BIE  vs.  other  numerical  methods  and  about  the  various  approximation 
methods  used  here.  Some  comments  regarding  extension  of  the  pre- 
sented methods  for  other  classes  of  boundary  value  problems  are  also 


given. 
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CHAPTER  TWO 

BIE  FORMULATION  AND  A BASIC  NUMERICAL  PROCEDURE 
2.  1 BIE  FORMULATION 

Consider  the  problem  of  determining  a function  u(x.)  throughout  a re- 
gion R,  i = 2 or  3 (for  two-dimensional  or  three-dimensional  regions,  re- 
spectively), whereu  satisfies  Laplace's  equation 

V2u  = 0 (x.)  € R (2.  1) 

and  R is  a simply  or  multiply  connected  domain  with  boundary  B.  The 
boundary  conditions  are  either  of  the  Dirichlet  type,  of  the  form 

u(x.)  = f(x.)  (x.)  £ B 

i i i 

or  the  Neumann  type,  of  the  form 

du(x.) 

= 6<x.)  (x.)  € B 

or  the  Churchill  type,  of  the  form 

u(x.)  = f(x.)  (x.)  £ B 

11  l 1 

du  (x.) 

^ * *<xi>  <*i»  * b2 

where  B - B,  + B_  and  — — is  the  derivative  in  the  direction  of  the  out- 

12  3n 

ward  normal  to  the  boundary.  Functions  f(x.)  and/or  g(x.)  are  prescribed 
on  the  boundary  B,  and  for  a well-posed  boundary  value  problem  f(x.)  and 
g(xj  are  not  simultaneously  prescribed  over  the  same  part  of  the  boundary. 


r 
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Now  for  any  two  sufficiently  smooth  and  non-singular  scalar  functions 
U and  V in  R,  it  is  well  known  from  Green's  theorem  [7]  that 

(U  v2V  - V v2U)dR  = C (U  & - V ^ )dB.  (2.2) 

jk  jr>  3 n on 

For  clarity,  we  now  consider  the  two-dimensional  and  three-dimensional 

cases  separately. 

A.  Two-Dimensional  Case: 

If  we  choose  U = u(t),  the  harmonic  function  we  are  trying  to  determine, 
and  V = log  r(t,w)  which  is  the  fundamental  singular  solution  to  Laplace's 
equation  in  two  dimensions  (r(t,w)  is  the  distance  between  any  two  points  in 
the  region),  we  will  find  that  Eq.  (2.2)  will  be  invalid,  because  V will  be 
singular  when  t coincides  with  w.  Therefore  we  construct  a small  circle 
B'  with  radius  r'  around  t as  shown  in  Fig.  2.  1, 
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and  now  Eq.  (2.2)  will  be  valid  for  the  altered  region.  In  addition,  because 
2 2 

both  V U and  v V are  always  zero,  Eq.  (2.2)  can  thus  be  written  as 

§B+B'*-U(q)  ~°^  ln'^  - lQg  r(t,q)  ] ds(q)  = 0 (2.3) 

1 

where  ds  is  the  differential  segment  of  boundary  B.  By  taking  the  limit  of 
the  integrals  over  B7  as  r7  -»0,  it  is  found  [7], 

lim  C , u(q)  al°S  * r/(t,q)d0(q)  = -2n  u(t)  (2.4) 

r7-0  J 

lim  lo?  r'(t-<!)  r/(t,q)d0(q)  = 0 (2.5) 

r7-0  9n 

Combining  Eqs.  (2.3)  (2.4)  (2.5)  gives  an  important  relation  which  we  will 
deal  with  later,  i.  e. , 

-*>  ■ -i  Sb  [»«)  ] - (*•« 

This  relation  expresses,  an  arbitrary  harmonic  function  in  terms  of  its 
boundary  values  and  boundary  values  of  its  normal  derivative. 

However,  as  mentioned  before,  we  know  that  for  a well-posed  boundary 
value  problem,  only  "half"  of  the  boundary  data  needed  in  Eq.  (2.6)  will  be 
given;  therefore  it  is  obvious  that  Eq.  (2.6)  is  not  a solution  to  a given 
problem  if  we  can  not  find  a way  to  get  the  other  half  of  the  boundary  data. 

In  fact,  how  to  find  the  other  half  of  the  boundary  data  is  the  basic  idea  of 
the  BIE  method. 

Let  p be  a point  lying  on  the  boundary  B,  then  by  drawing  a small  cir- 
cular bubble  around  p,  there  is  a new  boundary  B'  + B"  as  shown: 

I 

J 
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According  to  the  former  derivation,  we  know  that, 

■Wk  ■ -'<*<•>  <2'7» 

By  taking  the  limit  of  the  integrals  over  # and  B"  as  r'-O  it  is  found 

[7]  [8], 

lim  C u(q)  5l°S  I r'(p,q)de(q)  = (2tt-  a(p)Mp)  (2.8) 

, « on 

r'-O  J 

lim  C_M  log  r'(p,q)  r/(p,q)d0(q)  = 0 (2.9) 

r'-O  J 

where  a(p)  is  the  inner  angle  of  p.  Note  that  by  considering  a uniform  dis- 
tribution of  u over  the  region  mentioned,  <y(p)  can  be  expressed  as, 

«(P)  = 51°s^,~a)  ds<q>  <2-10> 

By  combining  Eqs.  (2.7)  (2.8)  and  (2.9),  we  arrive  at  the  basic  relation- 
ship of  the  BIE  method, 

u(p)  = ohp)  $B[u(q)  $1°S - lo8  r(P»9)  1 ds(9)  (2.1D 
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Eq.  (2.  11)  gives  us  a linear  relation  between  the  boundary  values  of 

du(p) 

the  harmonic  function  u(p)  and  those  of  its  normal  derivative  . 

The  essential  concept  of  the  BIE  method  is  therefore  the  following: 
Specify  the  known  boundary  data  in  Eq.  (2.  11)  and  solve,  by  one  of  the 
numerical  procedures  to  be  discussed,  for  that  part  of  the  boundary 
data  not  initially  prescribed.  Once  all  boundary  data  is  known,  u(t)  at 
any  point  t in  R can  be  obtained  by  means  of  Eq.  (2.6). 


B.  Three-Dimensional  Case: 

Let  U = u(t),  be  the  harmonic  function  to  be  determined,  and  V = 1/r 
(t, w)  which  is  the  fundamental  singular  solution  to  Laplace's  equation  in 
three  dimensions  (r(t,w)  is  the  distance  between  any  two  points  in  the 
region).  As  in  part  A,  Eq.  (2.2)  will  be  invalid  for  this  region,  so 
construct  a small  sphere  B'  with  radius  r'  around  t. 


Fig.  2.3 
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2 2 

Because  v U and  v V will  be  zero  for  the  region  R-R',  Eq.  (2.2)  can 
be  written  for  the  three-dimensional  case  as 

W"(q)  Tn  Tidr  - Tihr  ^ 1 d>«»  * 0 <2- 12> 

where  da  is  the  differential  area  of  boundary  B.  Taking  the  limit  of  the 
integral  over  B'  as  before,  it  is  found  that  [7] 


Um  $B,  u(q>  ^ 77^-  - r'2(t,q>d  r feu  <t) 

r -*0 

W bi>  ^ r'!(l',w  = 0 

Combining  Eqs.  (2.12)  (2.13)  and  (2.14),  we  thus  have 


(2.13) 


(2.14) 


■ ir  [ ±r>  - u(q)  fe  Tib 1 da,q)  <2- 15) 

Similarly,  let  p be  a point  on  the  boundary  B.  Construct  a small  bub- 
ble around  p thus  there  will  be  a new  boundary,  B'+B",  with 


/ y j'I  \ 

f 


Fig.  2.4 
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= *4n”  ^B'+B"  ^ " U^  ~%n  7(^T)lda(q)  <2' 161 


and 


lip  ^B„  u(q)  T/(p>q-)  r/2(p,q)dB  (q)  = («(p)  - 4^u(p) 

^ r'2,P'q|dB(q)  = 0 

■where  a(pV  the  inner  solid  angle  of  p,  can  be  expressed  as 


«<p>  = - $, 


da(q) 


>B  Qn  r (p,  q) 

Combining  Egs.  (2. 16)  (2.  17)  and  (2. 18)  we  thus  have 


(2.17) 


(2. 18) 


(2.19) 


u,p)  = $b[  Ti^T)  ^-  “<q)  t ^),da(q)  (2'201 

Now,  with  the  full  formulation  of  interior-and  boundary-boundary 
integral  equations  Eqs.  (2.6),  (2.11)  and  (2.15),  (2.20),  for  the  two- 
dimensional  and  three-dimensional  Laplacian  problems,  respectively, 
we  can  go  further  to  the  numerical  procedure. 


2.2  A BASIC  NUMERICAL  PROCEDURE: 

A Two-Dimensional  Case: 

The  boundary  B may  be  divided  into  N segments.  Thus  Eq.  (2.  11) 
can  be  written  as 


N 


Fig.  2.5 

Now,  there  will  be  N algebraic  equations,  each  corresponding  to 
the  boundary-boundary  integral  equation  (2.11)  for  point  p at  different 
segments: 


log  r(p, q)ds(q)  (2.21) 

J 

As  a first  approximation,  each  segment  may  be  modelled  as  a straight 

line,  and  the  function  u(q)  and  its  normal  derivative  u (q)  over  each 

n 

segment  may  be  considered  constant,  i.  e. , as  the  value  at  the  mid- 
point of  the  segment. 


a log  r(p. , q) 


Sbj  i^ds,q)tJI,uJ  Sb, 

“I  % Sb.  lo8r,pi-q)d’<q) 


alog  r(p.,q) 

— ds(q) 

3n 


ISisN 


(2.22) 


where  u and  u . are,  respectively,  are  approximate  values  of  u and 
J nJ 

its  normal  derivative  over  the  jth  segment.  For  conciseness,  we  can 
express  the  N linear  simultaneous  algebraic  equations  system  of  Eq. 
(2.22)  in  matrix  form: 


where 


[A]  {U}  = [B]  fUn} 

[A]  = [a..] 

[U]  = [u.  } 

[B]  = [by] 


[U  ) = [u  ] 
1 nJ  n. 

J 

Blog  r(p.,q) 


(2.23) 


lsiJSN 


ds(q) 


• t s 


a log  r(p.,q) 


k=l  K 


ds(q)  (2.24) 


b..  = ^ log  r(p.,q)ds(q) 


(2.25) 


Now,  it  is  apparent  that  we  can  construct  the  [A]  and  [B]  matrices 
numerically  from  Eqs.  (2.24),  (2.25)  (Appendix  1).  With  the  boun- 


dary data  that  is  given,  we  can  thus  solve  the  Laplacian  problem  as 


follows. 
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F or  the  Dirichlet  problem,  we  have 

(C)  = [B]  {Uj  (2.26) 

where  fC}  = [A]  fU} 

For  the  Neumann  problem,  we  have 

[A]  fU)  = {C}  where  fCt  = [B]  fUn}  (2.27) 

However,  note  that  the  Neumann  type  problems  do  not  have  unique  solu- 
tions, i.e.,  the  [A]  matrix  is  singular.  Therefore,  we  should  constrain 
arbitrarily  one  value  of  {U}  so  as  to  uniquely  obtain  all  the  other  values 
of  {U}.  An  easy  way  to  do  so  is  to  specify  one  u to  be  zero,  and  elimi- 
nate the  corresponding  row  and  column  in  [A]  and  the  corresponding  ele- 
ment of  fC}.  then  deal  with  the  reduced  system,  e.g. 

[A*]  {U*}={C*3  (2.28) 

where  [A*]  = [a.J  lSi,jsN-l 

{U*}  = {a.} 

{C*}=  [c.l 


Another  way  to  do  so  is  to  specify  arbitarily  the  value  of  one  element  of 
{U},  and  assume  the  corresponding  u to  be  unknown.  We  then  solve  the 


! 


! 


whole  system  as  a mixed  problem  as  described  below. 
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For  a mixed  problem,  we  have 


[A]  {U*}  = {C*} 


(2.29) 


where  [A*]  = [a*.] 

{U*}  = {u*3 

{C*}  = [B*]  fU*}  = [b*  ] Cu*.) 
‘ 1 n lj  nj 

For  a segment  where  u is  defined. 


lsi,  jsN 


a*.  = b..:  b * 

lJ  y*  y 


a..;  ur  = u u*. 
« J nJ  nj 


and  for  a segment  where  u is  defined, 

n 


a*  = a..  ; br.  = b..  : u*  = u.  ; u*.  = u . 

y y y y j j 7 nj  nj 


Eq.  (2.26)  Eq.  (2.28)  or  Eq.  (2.29)  can  be  solved  by  various  standard 
methods.  After  application  of  the  procedure  mentioned  above,  the  pre- 
viously unknown  boundary  data  is  determined. 

In  terms  of  the  discretized  boundary  and  the  entire  boundary  data, 
the  interior  field  equation  (2.6)  can  be  thus  written  as 


u(t)  = ~ [u.k  . - u k ] 
2tt  i li  n.  2i 

k^.  = log  r(t,q)ds(q) 


iSifiN 


(2.30) 


(2.31) 


(2.32) 


Mitt 


With  Eqs.  (2.30)  (2.31)  and  (2.32)  (Appendix  1),  u(t),  the  solution  of 


the  Laplacian  problem  concerned,  can  be  calculated  anywhere  of  in- 
terest in  the  region  R. 


B.  Three-Dimensional  Case: 

The  boundary  B may  be  divided  into  N segments.  Thus  Eq.  (2.20) 
can  be  written  as 


-U<P’I  k h T( £?)  da,q)+i  V T^q)  da,q) 


3 = 1 3 


3 = 1 J 


= y r SHiai  _J_  da(q)  (2 

L JlB.  r(p,q)  1 

j = l J 

As  a first  approximation  each  segment  may  be  modelled  as  a plane 


(2.3: 


triangle,  and  the  function  u(q)  and  it's  normal  derivative  u (q)  over 

n 

each  segment  may  be  considered  constant,  i.  e. , as  the  value  at  the 


centroid  of  the  segment. 


Fig.  2.6 
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Similar  to  the  two-dimensional  case,  there  will  then  be  N algebraic 
equations,  each  corresponding  to  the  boundary-boundary  integral  equa- 


tion (2.  19)  for  point  p at  different  segments: 

N N 

‘“"V  l Yn  r(iT7q>  da(q)  * I “j  $b.  f* 


da(q) 


N 


= l “nj  r,Pi'<!,da,q) 

J = 1 J 

In  matrix  form,  Eq.  (2.34)  may  be  written  as 
[A]  {U}  = [B]  [U  } 


where  [A]  =[  a„] 


fU}  = {u.} 
[B]  = [b..] 


IsisN 


[U  } = fu  } 
nJ  • n.J 
J 


and 


(2.34) 


(2.35) 


ij  J k £ da'ql  'hi  l k £ rl^T  da<q'  ,2‘ 36) 

J 1 k=l  K i 


b = C - 

ij  JBj  r(PL. 


q) 


da(q) 


(2.37) 


Once  the  [A]  and  [B]  matrices  have  been  calculated  numerically 
(Appendix  1)  with  half  of  the  boundary  data  prescribed,  the  other  half 
of  the  boundary  data  can  be  obtained  numerically  using  the  same  proce- 
dures as  in  the  two-dimensional  case  for  Dirichlet,  Neumann  or  mixed 


r 
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type  problems. 

With  the  entire  boundary  data  and  the  discretized  boundary,  the  in- 
terior-boundary integral  equation  (2.  15)  can  be  expressed  as 


u(t)  = 7-  [u  k - u.k  .] 
4tt  n.  2i  1 li 
1 


Cli  = Sb.  }n  r(t,q)  da(q) 


c = C -JL_ 

2i  JB.  r(t,q) 


da(q) 


IsisN  (2.38) 


(2.39) 

(2.40) 


As  in  the  two-dimensional  case,  u(t)  can  be  obtained  anywhere  of  in- 
terest in  the  region  R,  with  Eqs.  (2.39)  and  (2.40)  (Appendix  1). 


I 


CHAPTER  III 

IMPROVED  APPROXIMATION  METHODS 

In  the  preceding  chapter,  for  simplicity,  we  approximated  the 

function  u (q)  and  its  normal  derivative  u (q)  with  constants  over  each 

n 

segment.  Obviously,  in  many  cases  this  will  be  a very  crude  approxi- 
mation particularly  because  of  discontinuities.  Therefore,  improved 
approximation  methods  are  desirable.  There  are  various  classes  of 
approximating  functions  which  could  be  used,  for  example,  polynomials, 
trigonometric  functions,  exponentials,  and  rational  functions.  The  most 
widely  used  approximating  functions  for  well  understood  reasons  are  the 
polynomials  which  will  also  be  our  choice  in  this  work. 

As  a first  step,  we  will  discuss  what  polynomials  p(t)  are  good 
approximations  to  the  real  function  f(t),  i.e.,  u(q)  and  u^  (q)  in  our 
situation.  Before  that,  we  should  make  clear  what  is  meant  by  a "good 
approximation"  here. 

Definition: 

Let  f € Cn(a,b)  and  define  the  norm  ||  ||  of  f as 

II  i II,  = II £ lit  II f'  II +11  f"  ll+"-  + l|f(nlU 

where 

||f(j)||  = max  |f(j)(t)| 
t €(a,b) 

A function  g belonging  to  Cn(a,b)  is  a good  approximation  to  f,  provided 
|!  f - g jj  n £ €for  a sufficiently  small  €. 
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From  the  Weierstrass  approximation  theorem  [9],  which  says  for 
given  any  interval  (a,b),  any  real  number  € > 0,  and  any  real-valued 
function  f € Cn(a,b),  there  exists  a polynomial  p(t)  such  that  ||  f(t)  - 
p(t)  ||  L € for  all  t in  (a,b),  we  know  that  there  do  exist  polynomials, 
p(t),  which  can  be  good  approximations  to  the  real  function,  f(t).  The 
polynomial  we  are  interested  in  is  the  interpolating  polynomial. 

Let  a = t Z.t  Z.t  L . . .Lt  = b be  n + 1 distinct  points  on  a closed 
0 1 2 n 

interval  (a,b)  on  a real  axis  and  let  f (t)  be  a real  function  defined  at 

these  points.  It  is  obvious  that  there  exist  infinitely  many  polynomials 

p(t)  which  can  approximate  f(t)  and  interpolate  f(t)  at  the  points, 

t , t , t , t , . . . t , where 

U 1 c.  5 n 

p(t  ) = f(t  ) Osisn 


However,  there  is  only  one  polynomial  of  degree  £ n which  interpolates 
f(t)  at  these  n + 1 distinct  points.  In  terms  of  the  Lagrange  form,  it 
will  be 


where 


and 


p(t)  = f(tk)lk(t) 


Lk(t) 


n 

:in=0 

i^k 


(3.1) 


(3.2) 


W ■ 6*i 


0<j  5n 


1 (t)  defined  by  (3.2)  is  called  a Lagrange  polynomial  and  is  also  known 
k 


as  a "shape  function". 


Although  Lagrange  polynomials  are  very  easy  to  construct,  an  ob- 
jection to  this  type  of  interpolating  polynomial  is  that  when  the  number 
of  interpolating  points  becomes  large,  with  a consequent  high  degree  of 


Lagrange  polynomial,  calculation  and  evaluation  of  the  interpolating 
polynomials  become  costly  and,  due  to  the  propagation  of  errors,  un- 
reliable. Another  more  serious  objection  to  the  use  of  this  interpolat- 
ing polynomial  of  high  degree  is  the  fact  that  it  may  well  increase  the 
interpolating  error.  This  can  be  easily  seen  by  examining  the  error 
estimate  for  Lagrange  polynomial  p(t)  which  interpolates  f(t)  at  a = 

t L t L . . . L t = b.  It  is  found  [ 10]  that,  for  f € Cn+  *(a,  b) 

U 1 n 

r(n+l)  | 


II  £ - p II  * 


n+l 


4(n+l) 


for  f € C (a,b)  Osksn 


where 


f - p ||  4 1 + (^-)n]  gk(f,n) 


h = max  (t.  , - t.) 
0£i£n- 1 1+1  1 


h = min  (t.  , - t.) 
0*i«-l  1+1  ‘ 


gk(£,n)  = 


' 6»<£-  Tj4) 


3fb  - a)  i.  j, 
n " 


k=0 


k=l 


][f(k)H  litan 

(k  - 1)1  n 
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We  can  see  then  that  the  case  of  f € C (a,b),  if  (|  f ^ ^||  increases 


sharply  as  n becomes  large,  reducing  the  size  of  h may  prove  to  be  of 
help  even  if  the  knots  are  spaced  optimally,  e.g.  [ll] 


no 


(n+1) 


f - p 


(n+1) 


b - a 
2 (-T-) 


n+1 


where,  the  choice  of  Tchebycheff  points  for  the  interval  (a,b)  to  mini- 
mize the  approximating  error  are 

t.  = + a - (b  - a)cos[  2 (n++  \)TT^  0sj 

For  f € Ck(a,b)  (Osksn)  case,  [ 1 + (Zh/h*)11]  may  grow  too  fast  and 
outgrow  the  opposite  effect  due  to  (b  - a)/n  for  k = 0 and  ((b  - a)/n) 
for  lsksn. 

of  course,  error  estimates  are  sometimes  just  pessimistic;  how- 
ever, from  them  we  do  find  that  there  exists  a way  to  guarantee  the 
"smaller"  error,  even  for  small  n.  That  way  is  to  make  the  interval 
(a,  b)  small,  and  this  leads  to  the  "piecewise  interpolating  polynomial" 
idea.  Based  on  this  idea,  we  will  go  on  to  the  following  improved 
approximations . 


3.  1 SHAPE  FUNCTION  APPROXIMATION 
A.  Two-Dimensional  Case: 

As  mentioned  above,  in  order  to  have  convergence  as  the  number 
of  segments  increases,  it  would  be  necessary  to  use  piecewise  approxi- 
mating polynomials.  Within  each  subinterval,  the  shape  function  e.g., 
a Lagrange  polynomial,  can  be  constructed  linearly,  quadratically. 
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cubically,  etc.  The  choice  of  degree  depends  on  whether  we  want 
economic  advantage  or  other  advantages  which  can  be  achieved  by  in- 
creasing the  number  of  degrees  of  freedom.  Here  we  will  choose  the 
quadratic  shape  function  to  illustrate  and  apply  the  BIE  method. 

When  we  approximate  the  functions  u(q)  and  u^  (q)  by  shape  func- 
tions in  terms  of  Lagrange  polynomials,  i.e. 


u(q)  = u(qkUk(q) 


a„(q) * “„(qkuk(ql 

;r  each  segment,  where  for  the  quadratic  Lagrange  polynomial, 


ana  s 


ubstitute  into  the  two-dimensional  boundary-boundary  integral 


equation  (2.11),  it  becomes 

N k 


a(p)u(p)  =\  £ u(q^)  ^ M 


log  r(p,< 
Sn 


-ds(q) 


i=lj=l 


-I  I w L v,llog  rlp 

i=l  j=l  ^i 


,q)ds(q) 


However,  it  is  found  that  unless  the  boundary  segment  is  a straight 
line  or  we  are  going  to  model  the  curved  boundary  segments  as  straight 
lines  as  we  did  in  Chapter  2,  it  will  be  difficult  to  perform  the  approxi- 
mate integration  over  each  segment  if  we  do  not  know  the  geometric 
expression  for  the  boundary,  which  is  frequently  the  case.  In  turn, 
this  will  be  a drawback  if  we  want  to  use  a computer-aided  program 
to  make  the  geometric  partition  of  the  boundary  so  as  to  gain  efficiency. 
Therefore,  a parametric  transformation  should  be  used  initially  to 


A 
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transform  the  given  arbitrary  boundary  segment  into  a standard  bound- 
ary segment  to  which  a rule  of  approximation  integration,  like  the  Gaus- 
sian quadrature  we  shall  adopt,  can  be  readily  applied.  In  fact,  this 
can  be  done  by  using  the  so-called  "isoparametric  transformation"  [l2], 
where  we  use  the  shape  functions  used  to  approximate  the  functions  to 
also  approximate  the  boundary  geometry. 

For  this  purpose,  let  the  standard  boundary  segment  be  described 
by  the  curvilinear  coordinate  e as  shown  in  Fig.  3.  1.  The  quadratic 
shape  functions  related  to  this  coordinate  are 

1 , , 2 
N (e)  = - £e  + \e 

N2(e)  = 1 - e2  (3.3) 


3 . , 2 

N (e)  = fce  + £e 


The  x and  y coordinates  and  functions  u(q),  u (q)  can  thus  be  approxi 
' n 

mated  as 

I 
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x(e)  = Na(e)xa 


y(e)  = Na(e)ya 


(3.4) 


u(e)  = N0£(e)u°! 


, . ..a.  . a 

u (e)  = N (e)  u 
n n 


(3.5) 


& OL  CL  Qt 

where  x , y , u , u are  the  function  values  defined  at  the  nodes  of 
n 

the  boundary,  and  lsas3  for  the  quadratic  case.  It  is  also  found  that 


ds(q(e))  = J(e)de 


(3.6) 


where 


/_dx)  + 

' de  / ' de  J 


,/_dN^[e)  xa  2 , / dNtt(e)  a* 

\ de  J \ de  y J 


(3.7) 


Before  we  apply  these  shape  function  expressions  to  the  boundary- 
and  interior -boundary  integral  equations,  we  will  further  look  at  the 
boundary-boundary  integral  equations.  It  is  found  that,  by  using  Eq. 
(2.  11)  Eq.  (2.  10)  can  be  written  as 


^ [u(q)  - u(p)]  °S3nP,<:i  d5(q)  = ^ un(q)lo8  r(P»<l)ds(q)  (3.  8) 


Now  by  substituting  Eqs.  (3.5)  (3.6)  into  Eq.  (3.8)  for  each  specific 
point  p,  Eq.  (3.8)  can  be  formulated  as 


■ i riffiUtlHhiiniii  iUfi 
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N 3 &log  r(p  q(e)) 

I 1“  LN  (e) — s J,e,de 

J-l  «=1  Bj 


lsi,  j(a)  s2N 


Slog  r(p.,q(e)) 


J (e)de 


N 3 


= Z li  un3  a ^ Na(e)log  r(p^,  q(e))J (e)de 
j=l  OFl  n 


(3.9) 


where  N is  the  number  of  segments  and  B'  is  the  j standard  boundary 
segment  with  limits  (-1,  1).  In  short,  Eq,  (3.9)  can  be  expressed  as 


i ii 

j=1  a=i 


N 3 . , . 

y w «v. 

L L n ij 
j=i  ori 


(3.10) 


a ot 

Since  a , a'  b.  . can  be  calculated  numerically  (Appendix  2),  we 
ij  iJ  iJ 

can  thus  form  Eq.  (3.  10)  in  matrix  form  as 


[A]  {U}  = [B]  [Un} 


(3. 11) 


Thus,  for  either  Dirichlet  or  Neumann  or  mixed  type  problems  the  un- 
known half  of  the  boundary  data  can  be  obtained  by  the  procedure  de- 
scribed in  Sec.  2.2. 

With  all  these  boundary  data,  u and  u^  , lfims  2N,  the  two-di- 
mensional Laplacian  problem  can  thus  be  solved  by  the  use  of  the 
Interior-boundary  integral  equation  (2.6),  which  can  be  expressed  as 


d 
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“<*>  = 11  l i-1"-?, 

j=l  a=i 


k«  = C Nft(e)  51°g  de 

lj  Vj  3n 

2 j = ^ Na(e)log  r(t,q(e))J(e)de 


can  be  calculated  numerically  (Appendix  2). 


B.  Three-Dimensional  Case: 

Unlike  the  two-dimensional  case,  here  we  deal  with  the  integration 
over  curved  surface  boundaries  where  a boundary  segment  could  be  a 
curvilinear  triangle,  quadrilateral  ...  etc.  Here  the  boundary  seg- 
ment we  will  use  will  be  a quadrilateral.  Similar  to  the  two-dimen- 
sional case,  in  order  to  perform  the  numerical  integration  and  the 
automatic  mesh  generation,  the  parametric  transformation  is  used. 
Also,  we  let  B(  = (-1,  l)x(  -1,  1),  the  curvilinear  quadrilateral  with 
curvilinear  coordinates  e,  f,  be  the  standard  boundary  segment  as 
shown  in  Fig.  3.2. 


j 
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Based  upon  this  standard  boundary  segment,  we  could  easily  con- 
struct the  two-dimensional  shape  functions  by  forming  the  tensor  pro- 
duct of  the  appropriate  one-dimensional  Lagrangian  shape  functions 
which  we  used  before.  However,  due  to  the  larger  number  of  degrees 
of  freedom  involved  we  will  not  proceed  this  way.  Instead,  we  will 
choose  the  Serendipity's  elements  f 1 3]  for  our  usage,  where,  for  the 
biquadratic  approximation,  the  shape  functions  can  be  expressed  as 


N“(e,f)  = 1(1  + eQ)(l  + fQ)(e0  + fQ  - 1)  ls«s4 

Ntt(e,f)  = ^(1  + eQ)  (1  - f2)  a=  5,7 

Ntt(e,f)  = id  - e2)  (1  + f2)  a=6,8 


(3.15) 


where 


e 


0 


ee 


a 


f 


0 


ff 


a 
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For  the  isoparametric  transformation,  the  x and  y and  z geometric  co- 
ordinates and  functions  u(q),  u (q)  can  be  approximated  as 

n 

x(e,  f)  = N®(e,  f)x® 
y(e,f)  = N®(e,f)y® 

z(e,f)  = N®  (e,f)z®  (3.16) 

u(e,  f)  = N®(e,f)u® 

a (e.f)  = N®(e,f)u®  (3.17) 

n n 

where  x®,  y®  z®,  u®,  u®,  as  before,  are  the  function  values  defined  at 
the  nodes  of  the  boundary  nodes  shown  as  Fig.  3.2.  Also,  it  is  found 
that  (Appendix  2). 

da(q(e,f))  = J(e,f)de  df  (3.18) 


where 
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Similar  to  the  two-dimensional  case,  it  is  found  that  by  using  Eq. 


(2.  19),  Eq.  (2.20)  can  be  rewritten  as 

Og  !u'ql  - '“P)1  t 7 (^)da(q)  = ^ “n(q)  ri^qida(ql 


(3.20) 


By  substituting  Eqs.  (3.17),  (3.18)  into  Eq.  (3.20)  for  each  specific 


point  p,  Eq.  (3.20)  can  be  expressed  as 

V f uj(tt)  < Na(e  f)  ^ i 

L L*  i,#  N (e,i)  Sn  r(p.,q(e,f)) 


L V v 

j=i  a=i  B 


J(e,f)de  df 


- “">l>  I L t rlpTTqJeJ))  J<e’£)de  “ Ki.m*NN 


Ej 


l l <(a)  Sg,  NVf)  r(id(e.f))J<°' 
j = l a=l  “Bj 


f)de  df 


(3.21) 


where  N is  the  number  of  segments,  NN  is  the  total  number  of  nodes, 
and  B(  is  the  standard  boundary  segment  with  limits  (-1,  l)x(-l,  1). 
Eq.  (3.21)  can  be  expressed  in  the  abbreviated  form 


1 1 1 qi,ai  •?,  - q<n»qu  i ■ n ^ 

j=l  a=l  j=l  a=i 


(3.22) 


OL  OL 

After  a(  a.  b.  . have  been  calculated  numerically  (Appendix  2). 
ij  ij  ij 

Eq.  (3.22)  can  be  written  in  matrix  form  as 


[A]  {U}=  [B]  {Un} 


(3.23) 


which  can  be  solved  as  described  previously. 
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When  all  the  boundary  data  have  been  obtained,  the  three-dimen- 
sional Laplacian  problem  can  then  be  solved  by  use  of  the  interior- 
boundary  integral  equation  (2.  15) 


u(t) 


j=l  a=i 


n 2 j 


j(a).  a 

u k.  . 


(3.24) 


where 


= L N“<e-f|  dr  vd^J(e-£'ded£  <3-251 


B . 
J 


k?j  =SB,  N“le>£)  r(ti<e,£»J(e’£)ded£ 

j 

can  be  calculated  numerically  (Appendix  2). 


(3.26) 


3.2  SPLINE  FUNCTION  APPROXIMATION 

By  examining  the  shape  function  approximations  used  in  last  sec- 
tion, it  is  found  that  both  the  Lagrange  or  the  Serendipity  type,  guaran- 
tee continuity  of  only  the  function  itself  at  the  "knots",  where  a "knot" 
is  a point  at  which  the  piecewise- smooth  polynomials  are  joined  to- 
gether. Consider  the  one-dimensional,  piecewise-quadratic  shape 

function  of  Lagrange's  family,  for  example,  as  shown  in  Fig.  3.3, 

th. 

where  l.(t)  is  the  quadratic  Lagrange's  polynomial  over  the  i in- 


3 


terval. 


*6  *7 


% 


to 

-f<V) 


L’ct)  a -fitj 


Fig.  3.3 


*5 

<fcj  A3 


0 


It  has  been  noticed  that  there  exist  sharp  corners  at  these  knots 

(t  , t , t, , . . . ) due  to  the  lack  of  differentiability  of  l/(t)  at  the  knots. 

L *x  O 

Of  course,  if  the  approximated  function  is  smooth,  then  the  drawback 
of  lack  of  smoothness  will  cause  trouble  when  the  solutions  found  num- 
erically have  to  be  differentiated  at  the  knots.  Therefore,  it  may  be 
desirable  to  use  approximating  functions  which  are  not  only  continuous 
at  the  knots  but  also  have  continuous  derivatives  at  the  knots. 

There  are  various  ways  to  achieve  the  smoothness  requirements 
for  the  piecewise  polynomials.  One  way  is  to  use  values  of  derivatives 
as  well  as  the  function  itself.  However,  if  spline  functions  are  used, 
these  requirements  can  be  met  by  using  only  values  of  the  function  it- 
self [14],  [ 1 5 J , [l6],  [17].  As  Schoenberg  pointed  out,  the  fundamental 


I 


A 
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properties  of  spline  functions,  as  a class  of  piecewise  polynomial 
functions  with  continuity  properties  of  the  function  itself  and  its  deriva- 
tives, are  simple  and  the  mathematics  involved  elementary.  In  recent 
years,  there  have  been  numerous  works  with  regard  to  the  development 
of  the  theory  of  approximation  by  splines. 

Here,  we  will  look  at  the  cubic  B-spline  first.  It  is  known  as  the 
basic  spline,  because  Schoenberg  has  proved  that  every  cubic  spline 
can  be  written  as  a linear  combination  of  B-splines.  By  considering 

a partition  over  the  interval  a = t,  fit  fit  S ...  St  =b  with  evenly 

1 2 5 n 

spaced  subintervals,  the  B-spline,  B^(t),  is  defined  as 


(t  - t.  )- 

l-2 


t € (t.  ,,t.  ) 

1-2  i-I 


h3  + 3h2(t  - t ) + 3h(t  - t.^)2  - 3 (t  - t.^)3 


t € (t.  .,t.) 
i-I  l 


B.(t)  = ,<  h3  + 3h2(t.+  1 - t)  + 3h(t.+  J - t)2  - 3<t  ,+  1 . t)3 


(ti+2  - 


1 e <vw 


‘ e <‘i+r*i+2> 


V 0 


otherwise  (3.27) 


This  function  is  represented  graphically  in  Fig.  3.4. 


Fig.  3.4 


As  mentioned  before,  any  cubic  spline  s(t)  can  be  expressed  as  a 
linear  combination  of  cubic  B splines  (cf.  Fig.  3.  5): 

n+1 

s(t)  = ^ a.B^(t)  (3.28) 

i=0 

where  a and  a , the  two  additional  flexible  degrees  of  freedom,  are 
0 n+1 

caused  by  the  not  strictly  local  property  of  B splines. 
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To  determine  all  the  coefficients  a^,  besides  the  original  data,  e.g. 

s(t.),  lsi<n,  a choice  of  the  two  extra  degrees  of  freedom  should  be 

made,  e.g.,  s'(t  ) = f ' (t  ),  s'(t  ) = f'(t  ) or  s"(t  ) = f"(t  ) s"(t  ) = f"(t  ) 
linn  1 Inn 

etc.,  as  the  extreme  boundary  condition,  the  periodic  case  as  shown  in 
Fig.  3.6. 


It  is  found  that  Eq.  (3.28)  can  be  reduced  to 

n 

s(t)  = £ a-iBi(t) 

i=  1 1 

due  to  the  periodic  properties,  i.e.,  tQ  = t^,  t ^ = t^  and  BQ(t)  = B^t). 

Bn+l(t)  = Bl(t)- 

For  our  purpose,  we  will  use  the  periodic  spline  with  equal  intervals. 
Therefore,  in  terms  of  cubic  B splines  the  potential  function  u(q)  and  its 
normal  derivative  un(cl)  for  a specific  closed  boundary  can  be  approxi- 


mated as 
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N 

u(s)  = ) e.B(s  - s.) 

k ’ 


N 


u (s)  = Y f.B(s  - s.) 
n A.  J J 

J = 1 


(3.29) 

where  s is  the  curvilinear  coordinate  of  the  specific  closed  boundary: 


B (s  - s.) 
1 J 


B(s  - s.)  = < B (s  - s.)  > = < 


0 2h  (s  - s.) 

J 

(2  - (s  - s.)  )3 

L. 


(2  - 


hs(s  - s.)s2h 

3 

fs  - s.)  (s  - s.) 


- 4(1  - 


il  )3 


Os(s  - s.)sh 
J 


B(s.  - s)  1 l.  B(s.  - s)  s - s.sO 
^ J J J J 


(3. 30) 


and  h is  the  length  of  the  equal  intervals  between  nodes. 

Substituting  Eq.  (3.29)  into  the  two-dimensional  boundary-boundary 
integral  equation  (2.  11)  results  in 


N N 

tt(p  ) y e B(p  - s ) - C Y e.B(q(s)  - s.) 
jt*i  J 1 J JB  hi  3 3 


j = l 


91og  r(p.,q(s)) 
Sn 


ds 


N 


B(q(s)  - sjlog  r(p.,q(s))ds. 


(3.31) 
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After  careful  rearrangement,  it  is  found  that  Eq.  (3.31)  can  be  ex- 
pressed in  the  matrix  form 

[A]  {E}  = [B]  |T}  (3- 

■where 
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b.  . 
ij 


- sjlog  r(p.,q(s))ds 
(s  - sjlog  r(p.,q(s))ds 
( - s ) log  r(p^q(s))ds 


j - s ) log  r(p.,  q(s))ds 


aij,  . can  be  calculated  numerically  (Appendix  3). 
From  Eqs.  (3.24),  (3.30),  we  know  that 


where 


[D]  = 


ru]  = [D]  fE} 
{Un}  = [D]  CF} 


4 100 1 

1 400  0 

0 1 4 1 0 

0014 0 


(3.34) 


L.  1 000  14  j 

Therefore,  Eq.  (3.32)  can  be  rewritten  as 


J 


r 
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[A]  [D]"1  {U}  = [B]  [D]"1  {Un} 


i.  e. 


[A]  {U}  = [B]  {Un}  (3.35) 

Then,  by  the  similar  procedures  mentioned  in  Sec.  2.2,  all  the  neces- 
sary boundary  data  can  be  obtained.  It  should  be  pointed  out  that  a 
multiply  connected  region  can  also  be  handled  by  generalizing  the 
procedure  mentioned  above. 

After  all  the  u,  and  u have  been  determined,  the  two-dimen- 
sional  Laplacian  problem  can  be  solved  by  substituting  Eq.  (3.29)  into 
the  two-dimensional  interior -boundary  integral  equation  (2.6)  as 


e B(q(s)  - s.) 

J J 


Slog  r (t,  q(s))  ds 
dn 


N 

- Y f-B(q(s)  - s.)log  r (t,  q(s))ds] 

•>BjV  J 


(3.36) 


!•  6 • 


where 


(3.37) 
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C B,(.  S1°e  r(t.s(=))  de 

2 J 3™ 


+ r B,(. M°g  [(MW  a. 

V , 1 J a“ 

3+1 

r _ aiog  r(t,,fl(s))  ds 

JB  2 3 °n 

j-i 

+ ^ B (g  - S)  aiogr(t,q(S))  dg 
^B.  i J 9n 


k . = \ B (s  - s.)log  r(t,q(s))ds 
2J  JB  2 J 

3 


+ \ B (s  - s . ) log  r(t,  q(s))ds 
JB.  , 1 J 


3+1 


+ 


\ B,  (s.  - s ) log  r(t,  q(s))ds 

JB.  , 2 3 

3-1 


+ V B (s.  - s ) log  r(t,  q(s))ds 

JB.  , 3 

3-1 


can  be  calculated  numerically  (Appendix  3),  and  fe.} 


and 


can 


(3.38) 


(3.39) 

be 


obtained  via 


-43- 


3.3  TEST  PROBLEMS: 

3.3.1  Two-Dimensional  Case; 

(a)  Consider  the  temperature  distribution  along  the  boundary  of  a cir- 
cular disk  (see  Fig.  3.  7)  to  be  maintained  as  u(r  ,0)  = sin0. 


Fig.  3.7 

The  exact  analytical  solution  is  u(r,  0)  = (r/r  ) sin  0 and  it  is  desired  to 

o 

solve  for  the  boundary  heat  flux  by  using  the  BIE  method  with  piecewise- 
constant  approximation,  quadratic  shape  function  approximation  and  cubic 
spline  function  approximation.  Comparisons  are  listed  in  the  table  below. 


Approximation 

method 

Boundary 
points  used 

Average 
relative  error 

Computing 

time 

Constant 

10 

8.61% 

0. 25  sec. 

20 

2.  04% 

0. 64  sec. 

40 

0.49% 

0. 84  sec. 

80 

0.  11% 

3. 25  sec. 

Quadratic 
shape  function 

10 

0.  35% 

0.21  sec. 

20 

0.  10% 

1.18  sec. 

40 

0.  04% 

4. 1 3 sec. 

C ubic 

spline  function 

10 

l.  13% 

0. 94  sec. 

20 

0.11% 

3. 68  sec. 

40 

0.  01% 

3. 95  sec. 

All  the  relative  errors  listed  in  Sec.  3. 3 are  taken  as  the  absolute  value 
of  the  ratio  of  the  difference  between  the  analytical  solution  and  numerical 
solution  to  the  analytical  solution.  Regarding  the  boundary  solutions,  all 
the  relative  errors  at  each  boundary  point  are  then  averaged. 

Unless  otherwise  mentioned,  the  Gaussian  points  (Appendix  2)  used  in 
Sec.  3.3  are  4 and  4x4  for  the  quadratic  shape  function  approximation  and 
biquadratic  shape  function  approximation,  respectively. 

(b)  Consider  a circular  ring,  a multi-connected  region,  with  different 
constant  temperatures  along  the  inner  and  outer  boundaries,  as  shown  in 
F ig . 3.8. 


u,.«f 
d0  «o 
r\  *•  io 

rt 

F ig . 3.8 


The  exact  analytical  solution  is  u(r)  = u.  + (u  - u.)  In  (r/r.)/ln  (r  /r.). 

l o l 1 o 1 

We  solve  for  the  boundary  heat  flux  by  using  the  BIE  method  with  piecewise 
constant  approximation,  quadratic  shape  function  approximation  and  cubic 
spline  function  approximation.  The  comparisons  are  listed  in  the  table 
below: 
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Approximation 

method 

Boundary 
points  used 

Average 
relative  error 

Computing 

time 

Constant 

20 

2.  3% 

0. 6 sec. 

Quadratic 
hape  function 

20 

0.  1% 

1 . 2 sec. 

Cubic 

pline  function 

20 

1.6% 

1 . 4 sec. 

(c)  Consider  a rectangle,  kept  at  different  constant  temperatures  on  two 
opposite  surfaces  and  insulated  on  the  other  two  opposite  sides 

Y 


bcio 

CXm  ZO 

The  exact  analytical  solution  is  u = (u  + u )/2  + (u  - u )x/a. 

2 1 2 1 

It  is  desired  to  solve,  by  using  piecewise  constant  approximation  and  quad- 
ratic shape  function  approximation,  partly  for  the  boundary  temperatures  and 
partly  for  the  heat  flux  and  then  solve  for  some  interior  temperatures.  Re- 
sults are  shown  in  the  tables  below: 
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Piecewi.se  Constant  Approximation 
( 1 6 boundary  points  used) 

(0.43  sec.  computing  time) 


Interior 

Coordinates 

Analytical 

Solution 

Numerical 

Solution 

Relative 

Error 

6.0  0 

0. 2000 

0. 1965 

00 

0 0 

0. 5000 

0. 5000 

0 

-4.0  3.0 

0. 7000 

0. 7005 

0.  1% 

Quadratic  Shape  Function  Approximation 
(16  boundary  points  used) 

(1.06  sec.  computing  time) 


(a)  Consider  a cube,  kept  at  different  constant  temperatures  on  2 oppo- 
site surfaces  with  temperature  varied  lineary  along  the  other  four  adjacent 


surfaces,  according  to  u - u^  + (u^  - u^Jx/a  as  shown  in  Fig.  3.  10. 


We  solve  for  the  boundary  heat  flux  by  using  the  BIE  method  with  piecewise 


constant  approximation  and  biquadratic  shape  function  approximation.  The 
comparisons  are  listed  below. 


(b)  Consider  a plate  with  a hole  at  the  center  and  all  the  temperature 


Fig.  3.  11 


After  solving  for  the  boundary  heat  flux  by  using  BIE  method  with  biqua- 
quadratic  shape  function  approximation  for  different  Gaussian  points, 

4 and  5,  the  comparisons  are  listed  below: 


Quadratic  Shape  Function  Approximation 
(120  Boundary  points  used) 


(c)  Consider  a solid  sphere,  as  shown  in  Fig.  3.  12,  with  the  axisym- 
metrical  surface  temperature  distribution  u(r  , 0)  = cos  0 


Fig.  3.12 

The  exact  analytical  solution  is  u (r,  0)  = r/r  cos  0.  It  is  desired  to 

o 

solve,  by  using  quadratic  shape  function  approximation,  for  the  boundary 
heat  flux,  and  then  solve  for  some  interior  temperatures.  Results  are  giv 
en  in  the  following  table: 


CHAPTER  FOUR 


DISCUSSION  AND  CONCLUSION 

The  purpose  of  the  work  in  this  thesis  is  two-fold,  viz.  , 

(1)  to  apply  the  BIE  method  to  Laplacian  problems: 

(2)  to  supply  various  numerical  approximations  in  dealing  with  such 
applications. 

For  the  first  purpose,  it  is  obvious  that  the  BIE  method  has  a basic 
advantage  over  other  numerical  methods  like  the  FDM  and  FEM,  due  to 
the  reduction  of  dimension  by  one.  Thus  the  efforts  in  the  numerical 
modelling  and  preparing  the  geometrical  data  have  been  greatly  de- 
creased, especially  when  the  ratio  of  the  surface  area  to  the  volume  of 
the  domain  concerned  is  relatively  low.  However,  there  are  some  draw- 
backs with  the  BIE,  e.g., 

(1)  it  is  applicable  only  to  linear  problems; 

(2)  the  coefficient  matrix  involved,  i.e.,  [A]  and  [B]  shown  in 
Chapter  2 and  3,  is  a full  matrix,  therefore  not  like  the  banded 
and  sparse  matrix  involved  in  some  other  numerical  methods, 
like  FEM.  This  is  an  undesirable  computational  feature. 

However,  despite  these  drawbacks,  the  BIE  is  still  worthy  to  be  rec- 
ommended as  a useful  numerical  method  to  deal  with  boundary  value 
problems. 

As  to  the  second  purpose,  it  is  clear  that  one  can  not  say  in  advance 
which  numerical  approximation  is  superior  to  the  others.  The  term 
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"improved"  used  in  Chapter  3 simply  means  the  following: 

(1)  Instead  of  the  piecewise  step  function  approximation,  we  intro- 
duce the  higher  degree,  quadratic,  piecewise  interpolating  poly- 
nomial to  approximate  the  boundary  function.  From  the  error 
estimate  it  can  be  shown  that,  in  general,  the  latter  approxima- 
tion will  have  better  accuracy  than  the  former  with  the  same 
mesh.  Also,  by  using  the  isoparametric  transformation  as  ex- 
plained in  Chapter  3.1,  this  approximation  method  is  made 
more  useful. 

(2)  When  considerable  smoothness  of  the  representation  is  neces- 
sary or  desirable,  we  further  introduce  the  spline  function  ap- 
proximation. Here,  unlike  the  former  two  approximations,  con- 
tinuity of  the  derivatives  of  the  approximation  function  at  the 
knots  are  also  assured. 

Besides  these  two  facts,  there  are  various  other  factors  which  must  be 
taken  into  consideration  in  order  to  say  which  approximation  method  is 
improved  or  not.  Specifically,  consider  the  following: 

(1)  From  the  error  analysis  for  the  piecewise  quadratic  Lagrange's 
polynomial  approximation  [ll],  it  is  found  that 


where  f is  the  function  to  be  approximated,  p^  is  the  piecewise 
quadratic  Lagrange's  interpolating  polynomial,  and  h is  the  mesh 
length  interval  between  nodes.  An  error  analysis  for  the  cubic 


-53- 


spline  function  approximation  [2l]  leads  to 


II  * 


384 


where  s^  is  the  cubic  spline  function.  We  can  predict  that  the 
cubic  spline  function  approximation  will  have  the  faster  conver- 
gent speed  as  the  mesh  length  h decreases. 

(2)  In  order  to  achieve  the  same  accuracy,  the  quadratic  shape 
function  approximation  may  require  less  computing  time.  This 
has  been  shown  in  the  results  of  test  problems. 

(3)  The  effort  spent  on  the  programming  is  also  taken  into  consid- 
eration. The  piecewise  constant  approximation  is  relatively 
easy  to  deal  with.  In  contrast,  the  shape  function  and  spline 
function  approximations  are  comparably  complicated  to  program 

Therefore,  whether  a specific  numerical  approximation  is  "improved" 
or  not  to  any  analyst  depends  on  the  accuracy  needed,  computing  time 
available,  and  programming  efforts. 

In  Sec.  3.3,  we  have  posed  several  test  problem.  From  them, 
we  can  see  the  following: 

(1)  From  the  viewpoint  of  accuracy,  shape  function  approximation 
is  much  better  at  the  beginning,  while 

(2)  From  the  viewpoint  of  convergence  speed,  spline  function  ap- 
proximation is  faster  then  the  other  two  approximations; 

(3)  From  the  viewpoint  of  computing  time,  it  is  worth  while  to 
use  step  function  approximation  at  first  if  the  analyst  wants  a 
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rough  idea  about  the  solution  field. 

(4)  Accuracy  could  also  be  improved  by  increasing  the  number  of 
Gaussian  points  see  (Appendix  2). 

All  of  these  programs  were  run  on  an  IBM  370  Model  165  computer. 

There  are  further  considerations  with  the  different  approximations. 
When  we  use  step  function  approximation,  function  values  defined  are 
those  at  the  centroid  of  the  segment,  so  we  do  not  have  the  so-called 
"sharp-corner"  problem,  where  the  normal  derivative  of  the  potential 
function  is  discontinuous.  However,  in  the  shape  function  and  spline 
function  approximation,  we  do  have  the  sharp-corner  problem,  because 
in  these  two  approximations,  function  values  are  defined  at  the  boundary 
points  of  the  individual  segment.  To  deal  with  this  in  shape  function 
approximation,  it  should  be  remembered  that  only  one  unknown  can  be 
found  at  a nodal  point.  So  if  the  potential  function  is  an  unknown,  all 
the  different  values  of  its  normal  derivative  should  be  defined;  or  if 
potential  function  is  defined,  then  all  but  one  of  its  normal  derivatives 
still  needs  to  be  defined.  The  value  of  the  normal  derivative  can  be 
determined  via  the  boundary  geometry  and  the  potential  function  dis- 
tribution. Take  a two-dimensional  rectangle  as  an  example,  as 


shown  in  Fig.  4.  1, 


I 
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Fig.  4.  1 


' 


■where 

du2  dUj 

Un2  = dn  ~ dt  = a 
du  du_ 

u =-r±  = — = 0 

dn  dt 

and  t is  the  tangential  direction. 

As  to  the  cubic  spline  function  approximation,  because  of  its  in- 
herent smoothness  property,  the  function  to  be  approximated  should  be 
continuous  and  have  continuous  first  and  second  derivatives.  There- 
fore, in  dealing  with  the  sharp  corner  problem,  we  should  discard  the 
periodic  spline  approximation,  and  instead  we  should  use  the  ordinary 
spline  which  ends  at  these  points  where  discontinuity  occured. 

There  is  another  way  to  achieve  the  smoothness  requirement,  i.e.  , 
to  consider  the  piecewise  Hermitian  polynomial,  take  the  piecewise 
cubic  Hermites,  p(t),  for  example, 


N N 

P(t)  = v f(t  )h  (t)  + y f'(t.)h..(t) 

Lj  1 1U  i_,  1 ll 

i=0  i=0 

where  £(t)  is  the  function  to  be  approximated,  and  h.g(t),  h.^(t) 


(4.1) 


are  those  cubic  Hermitian  shape  functions  which,  similar  to  the  La- 
grangian  shape  function,  Eq.  (3.2),  have  the  following  properties: 


h (t.)  = 6.  . 
i0  j ij 


h,io(y  = ° 


h.  _(t.)  = 0 
J 


h;.  ,(t.)  = 6.  . 
il  J 1 1 

Thus  p(t)  has  the  continuity  of  the  first  derivative  of  f(t).  Substituting 
Eq.  (4.  1)  into  the  two-dimensional  boundary -boundary  integral  equa- 
tion and  interior -boundary  integral  equation,  we  can  solve  the  problem 
by  procedures  similar  to  those  described  in  Sec.  3.  l.A,  with  the  iso- 
parametric transformation.  However,  the  difficulties  are  two-fold. 


(1)  It  is  not  always  possible  to  obtain  accurate  values  for  the  deri- 
vatives of  the  known  function. 


(2)  Because  the  number  of  the  known  and  unknown  functions  has 
been  doubled,  the  size  of  the  resulted  coefficient  matrix,  [A] 
and  [B],  have  been  enlarged  from  N x N to  2N  x 2N,  where  N 


is  the  number  of  the  points  where  either  u or  u has  been  de- 

n 

f ined. 


Therefore,  we  did  not  consider  the  Hermite  Polynomial  method  in  this 


thesis,  but,  if  necessary,  it  can  be  implemented  readily. 

Finally,  as  is  pointed  out  in  the  beginning  of  this  chapter,  there  is 
no  absolute  advantage  of  any  numerical  approximation.  Therefore,  try- 
ing to  find  a "better"  approximation  in  order  to  improve  the  efficiency 
of  the  BIE  method  is  frequently  necessary.  Also,  it  is  possible  to 
generalize  the  methods  mentioned  in  this  thesis  to  the  boundary  value 
problems  governed  by  other  differential  equations,  e.g.,  Poisson's 
equation,  and  equations  which  govern  inhomogeneous  media  such  as 
bodies  with  inclusions  where  conditions  of  continuity  at  the  interface 


boundaries  need  be  considered. 


APPENDIX  I 


i 


CALCULATION  OF  ELEMENTS  OF  COEFFICIENT  MATRIX 
IN  PIECEWISE  STEP  FUNCTION  APPROXIMATION 

1.1  Two-Dimensional  Case: 


Fig.  A.  1.  1 
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lim  *V 
€'-0  Js' 


s°-€/  Blog  r (p. , q) 


.°-€' 


Bn 


• 9 T | 

ds(q)  = lim  \ — cos(0  - 

€'-0  V 1 


X)ds(q)  = 0 


s Blog  r(p  q)  a 

lim  \ ds (q)  = lim  \ n — cos(X  - 0 )ds(q)  = 0 

€.._0  -Js°+€"  r 2 


e"-o  +€ 1 


where  the  angles  in  the  expressions  above  are  as  shown  in  Fig.  A.  1.2. 
The  sum  of  these  two  integrals  always  has  a limit,  i.e.,  0,  no  matter 
how  the  €;  and  6"  tend  to  zero  independently  of  each  other.  Thus,  by 


the  definition  of  the  improper  integral,  we  know  that 


B log  r(p  , q) 

\ 7Z ds(q)  = lim  ( \ 

^B.  5 €'-0  ^ ^s" 


s°-€/  Blog  r(pi<  q) 


€ 

€M-0 


Bn 


ds  (q) 


S 

^s°+6" 


Blog  r (p  , q) 


Bn 


ds  (q)) 


= 0 


Therefore, 


a. . = 


N 

-5..  y 

ii 


0. 


1^1 

k^i 


ik 


= -rr 


1.1.2  Calculation  of  b. . : 
LL 


b..  = 


J 


log  r(p.,  q)ds(q) 


Again  there  are  two  different  conditions,  i ^ j and  i = j: 


(2.25) 


(i)  i 4 j 


b..  = 
n 


a o 

i = lil^  ^ log  r(p  q)ds(q)  + lim  { log  r(p., 
s-s°  ws  1 s-s°  1 


q)ds(q) 


Cr  rr 

= lim  \ log  r(p.,q)(-dr(q))  + lim  \ log  r(p.,  q)dr(q) 
0Jr"  i r-’0Jr  1 


= 2 lim  ^ log  r(p. , q)dr(q) 
r-0  ~r 


= 2 lim  [r(p.,q)log  r(p.,q)  - r(p.,q)]r 
r-0  1 1 1 r 


2[r+log  r+  - r+  - lim(r  log  r - r)] 


2 r+(log  r+  - 1) 


= As.[log(— 2^-)  - 1] 


cf . Fig.  A.  1 . 4. 


1.1.3.  Calculation  of  k. . and  k_ . 

li  2i: 


5log  r(t,q) 


ds(q) 


(2.31) 


;i  JL  1o8  r(t. 


q)ds(q) 


(2.32) 


Since  t is  not  on  the  boundary,  there  is  no  singular  case  and  we  can 


write  the  expression  for  k^  and  k,,.  directly  from  A.  1.  1.  1 (i)  and 


A.  1.1. 2(i),  respectively,  i.e., 


— — = - ~rr  (cos  u cos  a + cos  v cos  ft  + cos 

Sn  r(p.,q)  ^ 

= - "T  <5  ' a) 

r 

1 

= - — cos  6. 
r 

See  Fig.  A. 1.6. 

However,  it  has  been  noted  that 


Fig.  A.  1.7 


cos  0 


da'(q) 

da(q) 


and 

i _ de 

r2  “ da'(q) 

where  da'(q)  is  a segment  surface  of  the  sphere  with  radius  r, 

S 21 

tered  at  p.,  and  d0  is  the  solid  angle  confined  by  the  cone  p. 


w cos  y) 


cen- 

DEF,  See 


Fig.  A.  1. 7.  Thus 


da(q) 


a = L Gp.H 
1 

b = L Hp.I 
1 

c = L Ip.O 
1 

a+b+c 
5 = 2 


F ig.  A.  1. 8 


r 
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It  can  be  seen  from  Fig.  A.  1 . 8 that 

da"”  = L t 7^5) 


’ B* 

i 


^ — da(q) 

>B/  Bn  r(p.,q) 


and 


L & 7(^5  da,q)  * L • "2  ® ' •a)da(q)  ■ 0 


X5 

n 

i 


B*  r 

l 


That  is  to  say,  no  matter  how  £ approaches  zero,  the  integral  over 
this  plane  triangle  without  the  circular  area  with  radius  € and  center 
p.  will  always  have  the  same  limit,  0.  Thus  from  the  definition  of 
an  improper  integral  we  know  then 


da(q)  = lim  \ 


jb  dn  r(p.,q) 
i 


Therefore,  we  know  that 


e-o  -B*  Bn  r(pi'q) 


da(q)  = 0 


N 


sa 


a..  = - 6. . ) A0" 

11  11 L ik 

k^i 


=2t 


1.2.2  Calculation  of  b 


B r(p.,q) 


da(q) 


(2.37) 


(i)  i 4 i 


Fig.  A.  1.9 


In  spite  of  considerable  recent  activity  in  the  field  of  methods  for 
evaluating  multiple  integrals,  the  general  theory  is  relatively  undevel- 
oped. The  method  used  here  for  the  approximate  integrals  of  functions 
over  a triangle  (cf.  Fig.  A.  1.9)  is  the  64-point,  15th  degree  triangular, 
product  Gauss  formula  [ 1 8] . 

(ii)  i = j 

Taking  a right  triangle  into  consideration,  (Fig.  A.  1.  10)  we  know 


-71- 


Since  t is  not  on  the  boundary,  (cf.  Fig.  A.  1.  12)  there  is  no  singu- 
lar case,  and  we  can  write  the  expression  for  k directly  from 
A.  1.2.  1 (i). 


and  calculate  k_.  by  the  64-point, 

2i 

formula  as  in  A.  1.2.2(i). 


l 

15 


degree  triangular-product  Gauss 


APPENDIX  II 


CALCULATION  OF  ELEMENTS  OF  COEFFICIENT  MATRIX 
IN  QUADRATIC  SHAPE  FUNCTION  APPROXIMATION 


2.  1 Two-Dimensional  Case: 


2.1.1  Calculation  of  a?!  and  a'.: 

ii ii 


a. 

a. . 
U 


Slog  r(p  , q.(e)) 

N°(e)  . J J (e)de 

a n 

a log  r(p.,q.(e)) 

— — - — ^ J (e)de 


As  in  the  previous  Appendix,  we  will  have  two  cases  to  be  consid- 

th. 

ered.  (i)  p.  is  not  situated  within  the  j segment,  (neither  on  the  end 

th. 

point  nor  on  the  midpoint  of  the  j segment,  (cf.  Fig.  A.  2.  1). 


Fig.  A. 2.  1 


According  to  A.  1.  1.  1 (i),  we  know  that 


-73- 


Blog  rfp^q) 
3n 


d0 

ds 


Therefore,  a0*.  , a can  be  rewritten  as 


and 


1“  = C Ntt(e) 

U J_x 


d0..(e) 

— ±1—  de 
de 


*'•  = \ 

1J  J.i 


1 d0..(e) 

-ii—  de 


de 


0 ,(e)  = tan 


(y.(e)  - yfpJ) 
(x.(e)  - x(p .)) 


d0..(e) 

_U 

de 


y ! (e)(x.(e)  - x(p.))  - x'  (e)(y  (e)  - y(p.)) 
J - - J - - J ■■  - mmmm i 

2 2 
(x.(e)  - x(p.))  + (y.(e)  - y(p.)) 

3 1 3 1 


Let 


(e)  = Ntt(e) 


d0..(e) 

__u 

de 


G..(e) 


d0..(e) 

_U 

de 


we  can  thus  use  the  m-point  Gaussian  Quadrature  formula  [19]  to  cal- 
culate the  integrand  numerically  as 


a. . 
ij 


where 

(ii)  PL 


c**3 


Let 


= ) W„G..(e) 

L k ij  k 


w is  the  weight  function. 

K 

th. 

is  situated  within  j segment,  (cf.  Fig.  A.  2. 2) 


F ig.  A.  2 . 2 . 


a..  = a., 

ij  iJ 


lim  (Ntt(e)  - 

q.(e)-p. 


5..(a)a' 

ij  *3 


d0 . .(e) 

6..  (a)  — 4^ de 

ij  ae 


, d0  . . (e) 

...(a)' 

ij  de 


(N  (e)  - 6..(o)  )(y'  (e)(x.(e)  - x(p.))  - x'  (e)(y.(e)  - y(p.))) 

n J 1 i J J i 


(x.(e)  - x(p.))2  + (y.(e)  - y(p.))2 
i i ) 


which  is  an  indeterminate  form.  However,  by  using  l'Hopital's  rule 
twice,  we  found  that 


1 d0..(e) 

iim  (N  (e)  - 5..(o)  — ^ — 

/x  lj  de 

q.(e)-p. 


(N  (e))"(y/.  (e)(x.(e)  - x(p.))  - x'  (e)(y.(e)  - y(p.))) 


2(Na(e))'  (y''(e)(x.(e)  - x(p.))  - x''(e)(y.(e)  - y(p.))) 


(N  (e)  - 6i.(a)(y'.'(e)x' (e)  - x';(e)y'(e)) 


D = 2((x^(e)  - xCpJJxj'  (e)  + (x'  (e))2  + (y.(e)  - y(p  ))y"(e)  + (y'  (e))2) 

Therefore,  we  know  in  fact  that  there  is  no  singularity,  so  that  we 


can  calculate  a . as  we  did  in  A.  2,  1 . 1 (j);  i.  e. , let 
ij 

tt  x d0  (e) 

F..  (e)  = (Na(e)  - 5. .(a)'  

ij  ij  de 


and 
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- a 

a. . 


m 


c=l 


. a 

wF(e) 
k ij  k 


2.  1.2  Calculation  of  bQ. : 

il 

1 

b“  = C Ntt(e)log  r(p.,q.(e))J(e)de 
iJ  J_i  i J 

As  before,  we  have  two  different  conditions  to  be  handled. 

(i)  p.  is  not  situated  within  j*  segment  as  shown  in  Fig.  A.  2.  1.  There 
is  no  possibility  to  have  a singularity;  therefore  let 


H°i  (e)  = Na(e)log  r(p.,q.(e))J(e) 

ij  1 3 


and  using  the  m-point  Gaussian  quadrature  formula,  we  have 

m 

b“  = ) w H*  (1  ) 

ij  L k ij  k 

k=  1 
th 

(ii)  p.  is  situated  within  j segment  as  shown  in  Fig.  A.  2.  2: 


i 

btt.  = C Ntt(e)log  r(p.,q.(e))J(e)de 
1j  J_  i 1 J 


We  have  three  cases  to  discuss: 


(a)  Let  j(  1)  = i,  as  shown  in  Fig.  A.  2.  2 (a) 

0 r 1 

b..  s C Ntt(e)log  r(p.,q.(e))J(e)de  + \ Ntt(e)log  r (p. , q.  (e))J(e)de 
1J  J 1 J «Jq  1 J 

The  second  part  of  the  right-hand  side  is  free  of  singularity,  so  we 
treat  it  as  in  A.  2.  1. 2 (i);  i.  e.  , let 


d 
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H“  (e)  = Ntt(e)log  r(p.,q.(e))J(e) 


As  to  the  first  part,  we  let 


e = f - 1 


thus,  de  = 2 f df 
0 


V Ntt(e)log  r (p. , q.  (e))J (e)de  = C N°i(e(f))log  r (p. , q.  (e(f))J (e(f)) ( 
^ 1 J i 1 J 


-2f)df 


where 


lim  log  r(p.,q.(e(f)))  • df  = lim  log  f * f = 0 


i J 


f-0 


while  N°*(e(f))  and  J(e(f))  are  finite  within  the  range  f £ (-1,0);  there- 
fore we  can  let 


_a 


H..(f)  = N (e(f))log  r(p.,q.(e(f)))J(e(f))(-2f) 
ij  1 J 


and  use  the  m-point  Gaussian  quadrature  formula 


m a. 


= V w,(H..  (f  ) + H“  (e  )) 

1J  L,  k 1J  k lj  k 

k=  1 


where  f € (-1,0)  and  e € (0,1). 


Similarly,  we  can  use  the  following  expression  for  the  other  two 


cases : 
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(b)  j(2)  = i,  as  shown  in  Fig.  A.  2. 2 (b). 
Let 


-e  = f 
2 

e = g 
therefore 


for  e € (-1,0) 
for  e €(0,1) 


and  let 


u 

\ Na(e(f))log  r(p.,q.(e(f)))J(e(f))(-2f)df 

1 

+ ^ Na(e(g))log  r (p. , q.  (e(g)))J(e  (g))(2g)dg 

J0  1 J 


H“(f)  = Na(e(f))log  r(p.,q.(e(f)))J(e(f))(-2f) 


K W - - ** 


Thus 


m 

b“  = t W (Ha.  (f  ) + h“  (g  )) 
lj  Li  k ij  k lj  k 


k=  1 


where  f € (-1,0)  and  g €(0,1). 


(c)  j (3 ) = i,  as  shown  in  Fig.  A. 2. 2 (c). 
Let 


-e  = 


2 

g 


1 


for  e €(0,1) 


therefore 
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0 

b..  = C Na(e)log  r(p.,  q.(e))J(e)de 
U J_j  i J 

+ C Ntt(e(g))log  r(p.,q.(e(g))J(e(g))(2g)dg 

J0  1 J 


and  let 


H®  = Na(e)log  r(p.,q  (e))J(e) 


H.“  = Na(e(g))log  r(p.,qj(e(g)))J(e(g))(2g) 


Thus 


m 

b“  = f W,  (H..  (e,  ) + H®  (g,  )) 
1J  k=l 


ij 


ij 


where  e € (-1,  0)  and  g € (0,  1). 


2.  1.3  Calculation  of  k®.  and  k®. : 

Ll ll- 


\ N®( 

J-1 


3 log  r(t,q.(e)) 
e)  — ^ J (e)de 


on 


1 

k2.  = ^ N®(e)log  r(t,q(e))J(e)de 


(3.  13) 


(3.14) 


d0.(e)  y'(e)(x.(e)  - x(t))  - x'  (e)(y.(e)  - y(t)) 

^ — J . , . — J-  — - _ . J J 

6 (x  (e)  - x(t))2  + (y.(e)  - y(t))2 

J J 

2 . 2 Three-Dimensional  Case: 

2.2.  1 Derivation  of  the  parametric  representation  of  boundary  sur- 


face; Since  a surface  can  be  represented  in  parametric  form 


r(e,f)  = x(e,f)i  + y(e,f)j  + z(e,f)k 


where  _r(e,f)  is  a vector  function  of  two  independent  real  variables,  e 
and  f,  the  so-called  parameters  or  curvilinear  coordinates,  and  is 
single-valued  and  continuous  in  a simply-connected  bounded  region  of 
the  e-f  plane.  It  is  found  [ Z 0]  that 


da(q(e,f))  = J(e,f)de  df 


where 


J(e,f)  = V EG  - F 


E = r e • re 


F = re  * rf 


G = rf  • rf 


Thus 
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2 2 2 
J(e,f)  = V (yezf  - yfze)  + (xezf  - x^)  + (x^  - x^ 


*2  -'.-2  *2 

N!  + ^ + N3 


The  unit  normal  vector  of  surface  concerned,  n can  be  expressed 


r x r , 
— e —f 


r x Lc 
—e  i 


IN  | JN|  IN  | 

J (e,  f)  ~ + J (e,  f)  + J(e,  f)  - 


= ni  + nj  + nk 
x—  -Y-  z— 


2.2.2  Calculation  of  a.,  and  a'.: 

*J  ij 


1 1 


“ = [ C Ntt(e,f)  ~ — 1 ■— J(e,f)de 

ij  Sn  r(p.,q  (e.f)) 


1 1 


[ ~~  ^—7 — 77,  J(e,f)de 

1J  J-i  3n 


The  two  different  cases  we  considered  are,  (i)  p^  is  not  situated 


within  the  j segment,  as  in  Fig.  A.  2. 4. 


f 


r(p.,q.(e,f)) 
i J 


Sn 


2 

r 


(p.,q.(e,f)) 
1 J 


1 

r(p.,q.(e,f)) 
i J 


[(x.(q(e, 

J 


+ (Yj(q(e,f))  - y(p.))n^ 


F“  (e,  f)  = Na(e,f) 


_a_ 

an 


1 

r(p.,q.(e,f)) 


G..(e,f) 


_d_  1 

9n  r(p. , q. (e, f )) 

* l t 


J(e.f) 


- x(p.))n 
1 ■* 


(zj  (g(e,  f)) 


J(e,f) 


-84- 


We  can  thus  use  (ml  x m2) -point  Gaussian  quadratic  formula  to  calcu- 
late 

a , 

a. a. . as 
ij  U 


ml  m2 

k= 1 1=1 


ml  m2 

i'ij  -l  wk  z 

k=l  1=1 


lilt- 

Z wiGij(vfi) 


(ii)  p.  is  situated  within  the  j segment,  i.e.,  to  be  one  of  the 
nodes  of  the  segment,  either  corner  points,  or  midpoints  as  shown 
in  Fig.  A.  2.  5. 
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1 1 


L®  = C { (N°(e,£)  - 6.  ) [ 

ij  J.!  lJ(«) 


Dx.  .N*  + Dy.  ,NJ  + Dz.  .N* 

— U-J ^-3  ] de  df 


r (p.,q.(e,f)) 
i 3 


where 


Dx.  . = x.(q(e, f))  - x(p.) 
ij  3 1 


Dy_  = y^(q(e,f))  - y(pj 


Dz.  . = z. (q(e , f))  - z(p.) 
ij  3 i 


The  two  typical  cases  are  shown  in  Fig.  A.  2.  5 (a)  and  (b).  We  examine 


only  these  two  cases;  all  the  other  cases  can  be  handled  similarly. 


(a)  F or  region  I of  Fig.  A.  2.  5 (a),  we  can  consider  the  polar 


coordinates  as  in  Fig.  A.  2.6. 


Fig.  A.  2. 6 


Thus 


A 
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_ 2sec0  (Ntt(e,f)  - 5.  ,(a)(Dx.  ,N*  + Dy  .N* 

(a“)  =-  C C iJ-i LLJ. 

lJ  J0  4)  r (D,.q.(e.f) 


r (p.,q.(e,f) 
1 J 


+Dz.  .N* 
1 1 3 


dr  d0 


where 


and  let 


e = r cos0  - 1 


f = rsin0  - 1 


F.°i  (r.8) 


-(N  (e,f)  - 5.  .(a))(Dx.  .N*  + Dy.  .N*  + Dz.  .N*) 
- iJ  *J  1 iJ  2 ij  3' 


It  is  found  that 


*1*  (r’6)  = “f 

r -*0 

which  is  indeterminate.  However,  by  using  l'Hopitai's  rule  twice  we 


lim  F . (r , 0)  = lim 


-[  {Na(e,f)  - 6.  ,(a))(Dx.  .N*  + Dy.  .N*  + Dz  .N*]" 
L] LJ  1 i]  2 li  3 


where  the  limit  of  the  numerator  part  will  be  some  finite  number. 
Therefore,  we  know  that  there  is  no  singularity,  thus 


ml  m2 

(ay>i=I  l V?j  <V6i> 

k=l  1=1 
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where  r 6 (O,2sec0)  and  0 € (o,^tt).  Similarly,  we  can  conclude  that 


(a.  ) 


ij  II 


ml  m2 

w 


k=l  1=1 

where  r £ (O,2csc0)  and  0 £ so 


l "kl 


-ot  /-a  v , ,-a  v 

a.  . = (a,  .)  + a.  . 

ij  i j I ij  II 


(b)  Similarly,  by  considering  polar  coordinates,  for  region  I of  Fig. 
A.  2.  5 (b),  we  have 


■1. 


tan  (2)  sec0 

(a“)  = - C C F.a(r,0)drd0 

lJ  1 Jn  Jn 


U 


where  F ^ (r,0  ) is  defined  as  in  (a),  but  here  we  have 


e = rcos© 


f = rsin0  - 1 


Similar  to  the  procedure  in  (a),  we  can  have 


ml  m2 


”kl  "i  <vv 


ij  I L,  k 

k=l  1=1 


1, 


where  r € (0,  sec0)  and  0 £ (0,tan  (2)). 


“01  “0 1 

Also,  we  car,  have  (a  .)  and  (a.  in  the  same  form,  except 

ij  II  ij  III 


-1 


r € (O,2csc0),  0 € (tan  (2),rr-  tan  (2))  for  region  I 


-1, 


r € (0,  -sec0 ),  0 € (tt  - tan  (2),rr) 


for  region  II 


then  we  can  get 


J 
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a“.  = (a%  + (a®)  + (a“.) 

ij  ij  I ij  II  ij  III 


2. 2.  3 Calculation  of  b. 

LL_ 

I 1 


b“.  = C C Na(e,f)  l— 

ij  Xi  J.i  r(p.,q.(e. 


f)) 


J(e)  de  df 


The  two  different  cases  will  be 


.th 


(i)  p.  is  not  situated  within  j segment,  as  shown  in  Fig.  A.  2.  4.  Let 


1 


H“  (e.f)  = N“,e,£)  r(p_>qi(e>fl) 


J(e) 


1 J 


and  we  thus  have 


a 


ml 


m2 


w. 


ij 


c=i  r=i 

.th 


1 *iH“(ek’£i 


(ii)  p.  is  situated  within  j segment,  as  shown  in  Fig.  A.  2.5.  As  in 
2.2.2,  we  just  consider  the  typical  case.  Fig.  A.  5(a)  and  (b). 

(a)  By  considering  the  polar  coordinates  for  region  I as  shown  in 
Fig.  A.  2. 6 we  have 

(b®>  = ^ [ N“(e,f)J(e,f)  dr  d0 

ij  I o0  J>0 


where 


e = rcosG  - 1 


f = rsinG  - 1 


Let 
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H®  (r,  0)  = Na(e,f)J(e,f) 

CL 

It  is  obvious  that  H . is  free  of  singularity,  therefore 
U 


ml  m2 


(ba)  = y w,  y w H®. (r  0 ) 
ij  L k L 1 IJ  k 1 

k=l  !=! 


where  r £ (O,2sec0)  and  0 0 (0,^tt). 

q tt 

Similarly,  (b.  .)  has  the  same  form  as  (b.  .)  . except  that  r €(O,2csa0), 
ij  II  ij  I 

0 € (£tt,^tt)»  and  thus 

.a  a . ,,  a , 

b.  . = (b.  .)  + (b.  . 

ij  iJ  I ij  II 

(b)  Similar  to  the  procedure  described  in  (a),  we  have 

b“  = (b“.)  + (b®)  + (b“) 

ij  ij  I ij  II  ij  III 

OL  CL  CL 

where  (b.  (b  .)  and  (b  .)  have  the  same  form  as  shewn  in  (a), 
ij  I ij  II  ij  III 

except  that 

e = rcos0 


and 


f = rsin0  - 1 


r £ (0,  sec0),  0 0 (0,  tan  *(2))  for  region  I 


r £ (O,2csc0),  0 £ (tan  *(2),tt-  tan  *(2))  for  region  II 


r £ (O,-sec0),  0 € (tt  - tan  *(2),rr) 


for  region  III 


2.2.4  Calculation  of  k, . and  k„.: 
li  2j 

1 1 


“ij  3 5.1N“(6-£)  * r(t.  q^(e,  f,)  ** 


1 1 


*2j  = S.j  S_!  N (6,f)  r(t,q.(e,f))  J(e’f) 


Fig.  A. 2.1 


Since  t is  an  interior  point,  (cf.  Fig.  A. 2.1)  there  is  no  singular 

a o. 

case.  We  can  thus  calculate  k , k_.,  by  using  (ml  x m2)  -point  Gaus 

* J ^ J 

sian  quadrature  formula,  as  in  t lose  procedures  described  in  A.  2.2.2 
(i)  and  A. 2. 2.  3 (i);  i.e. 


ml 

m2 

■I 

wk  z 

k=l 

1=1 

r1 

m2 

= I 

k=  1 

1 = 1 

where 
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F?(e,f)  =Na(e,£)  ^ \ ...  J(e,f) 

j dn  r(t,  q.(e,  £)) 


Na(e, f)  ne.i) 


J 


and 


9n  r(t,q^(e,f) 


J(e,f)  = - 


1 


3 “ [(x.(q(e,f))  - x(t))n  + 

r (t,  q.(e,  £))  J 


(yj(q(e,f))  - y(t))n^  + (zj(q(e,f))  - z(t))nz] 


1 


APPENDIX  III 


CALCULATION  OF  ELEMENTS  OF  COEFFICIENT  MATRIX 
IN  CUBIC  SPLINE  FUNCTION  APPROXIMATION 


3.  1 Calculation  of  a..: 

U 


a..  = 
U 


3n 


ds 


9log  r(p.,q(s)) 


bi(s  ■ si)  ' - 

JBi+l  1 3 3 n 


B.(S.  - .)  5l°g  ''Pi-*'""  . 

2 1 ^ dS 


-Cr  Bis.-s,  al°«r(lVq(s|) 
JBi-2  1 J S d9 


(3.33) 


There  are  various  ways  to  calculate  these  integrals  numerically.  One 
convenient  way  for  the  general  case  which  we  will  illustrate  is  to  adopt  the 
parametric  transformation  of  geometrical  coordinates,  i.e.,  to  approxi- 
mate each  of  the  two  groups,  (B.,B.  ) and  (B.  _,B.  ) in  the  expression 

jj+l  J“* 

of  a..,  as  we  did  in  Chapter  3.  l.A.  (Consequently  we  should  divide  each 
closed  boundary  into  an  even  number  of  intervals  to  carry  out  this  kind  of 
approximation.  See  Figures  A.  3.1  and  A.  3.2.) 


Fig.  A.  3.  1 


Taking  group  I(B.,B^^)  as  an  example, 

3 0 d9..(e) 

'Vr-^V'^-V 

l d6..(e) 

-$  BL(s<e)  - 3j)  — ^ de 


■where 


Bjfsfe)  - sj  = (2  - 


B (s(e)  - s.)  = (2  - 
^ J 


(s(e)  - s ) 3 


(s(e)  - s.)  , (s(e)  - s ) 3 

r 4(1  - r !-) 


log  (2 


c(c  - b + a)  - 2c  + b)  I 
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a-*!>3  1-2  - 3 12 


(x  - x 


r + (y  - y V ] 

b = (x1  + x3  - 2x2)  (x3  - x1)  + (y1  + y3  - 2y2)  (y3  - y1) 


/ 1 3 2 2 1 3 2 2 

C = (x  + x - 2x  ) + (y  + y - 2y  ) 


Then  we  can  calculate  (a  . )j  in  the  manner  shown  in  A.  2.  1.  1 for  both  the 
nonsingular  and  singular  cases.  Similarly,  we  can  express  (a  ) in  the 

l ij  II 

same  way,  except  where  s.  - s(e)  = \ J(e)de,  and 

J Je 


a..  = (a..)  + (a..) 

U ij  I ij  II 


3.2  Calculation  of  b : 

11_ 


b-j  = - B2(s  - s.)logr(p.,q(s))ds 


' B (s  " s.)logr(p  ,q(s))ds 

1 J 1 


j+1 


■ B_(s  - s)logr(p  ,q(s))ds 

JSj.l  2 j 

“ - s)logr(p.,q(s))ds 


(3.34) 


Still  taking  two  groups  (B.,Bj+J)  and  (B^  2,Bj  as  in  A-3*  I and  Fig 

A.  3.  1 into  consideration,  we  thus  have 
0 

(b.Jj  = - ^_i  B2(s(e)  - s.)logr(p..q(e))J(e)de 


A 

C B (s(e)  - s.)logr(p.,q(e)J(e)de 
Jn  1 J 1 


1 


where  s(e)  - s.  is  defined  in  A.  3.1. 


Fig.  A.  3.2 

Without  consideration  of  singular  cases,  k and  k_.  can  be  calculated  by 

‘■3 

the  same  procedures  described  in  A.  3.  1 and  A.  3.2  respectively.  By  using 
a procedure  similar  to  that  described  in  A. 2. 1.2,  we  can  calculate  (b^)j 
for  both  the  nonsingular  and  singular  cases.  Then  by  combination  as  be- 
fore, we  have 


b„  - flO,  + <b..)IT 
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3.3  Calculation  of  k, . and  k_.: 


V $b.b2  «» 


+ Sb  b ai°g  de 

J+1  1 J an 


B2(Sj  - s) 


log  r(t,< 


+ L B (s . - s)  3-lQg 
-Pj_2  1 J dn 

k2j  = ^B.  B2(s  " Sj)log  r(t»9(s))ds 

+ \r  Bl(s  ' s;)lo6  r(t,q(8))ds 
'H+l  J 


+ Ki 


- s ) log  r(t,q(s))ds 


+ Bl(si  ‘ s)1°S  r(t* 9(s))ds 

J j-2  J 


(3.38) 


(3.39) 
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